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Abstract 

The peak of the two-particle Bose-Einstein correlation functions has a very inter- 
esting structure. It is often believed to have a multivariate Gaussian form. We show 
here that for the class of stable distributions, characterized by the index of stability 
< a < 2, the peak has a stretched exponential shape. The Gaussian form corre- 
sponds then to the special case of a = 2. We give examples for the Bose-Einstein 
correlation functions for univariate as well as multivariate stable distributions, and 
check the model against two-particle correlation data. 

Key words: correlations, elementary particle physics, heavy ion physics, statistical 
analysis, Levy-stable distributions 



1 Introduction 



The structure of the peak of two-particle Bose-Einstein correlation functions 
(BEC-s) is a very interesting physically measurable quantity. It carries infor- 
mation about the fraction of long-lived resonances [1], hence signals possible 
Ua{^-) symmetry restoration in high energy heavy ion reactions [2], it can be 
utilized to determine a generalized Hubble constant for the transverse flow in 
high energy heavy ion physics [3] and it is also sensitive to a possible par- 
tial coherence in the source [4], for other examples see the recent reviews of 
refs. [5,6]. 
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Originally, the method of Bose-Einstein correlations was invented by the radio 
astronomers R. Hanbury Brown and R. Q. Twiss (HBT), who applied it to 
determine the angular diameters of main sequence stars [7,8]. They referred 
to the method as intensity interferometry. In particle physics, an analogous 
effect was discovered independently [9] by G. Goldhaber, S. Goldhaber, W. - 
Y. Lee and A. Pais (GGLP), who observed an enhancement of identical pions 
with small opening angles in p + p reactions. Essentially, the enhancement of 
identical bosons with small momentum difference is due to their Bose-Einstein 
(BE) statistics. The effect is frequently referred to as either the HBT, or the 
GGLP effect, or intensity interferometry, or simply as BE correlations (BEC). 

The detailed shape analysis of the BEC is rather complicated but impor- 
tant, because the shape carries information about the space-time structure 
of the particle emission process. Recently, two of us proposed [10] a model- 
independent method for the analysis of the shape of this function, based on 
an experimentally preferred weight function and a complete orthogonal set of 
polynomials, where orthogonality is defined with respect to a preferred weight 
function. In the case of approximately Gaussian correlation functions, the 
method results in the Edgeworth expansions, and the complete orthogonal set 
of functions are the Hermite polynomials. In the case of an approximately ex- 
ponential shape, the expansion is given in terms of the Laguerre polynomials. 
For approximately spherical distributions the expansion can be formulated in 
terms of spherical harmonics (see ref. [10] for further details). 

The Edgeworth and the Laguerre expansions are very general methods to 
characterize the two-particle correlation functions. They depend only on the 
following two 

experimental conditions: 

i) The correlation function tends to a constant for large values of the relative 
momentum Q. 

ii) The correlation function deviates from its asymptotic, large Q value in a 
certain domain of its argument. 

For simplicity, this domain, where the correlation function deviates from its 
asymptotic value, with other words, the location of its non-trivial structure, 
is assumed to be close to Q — 0. 

This implies that the Edgeworth, Laguerre (and similar) expansions can be 
applied to all kinds of data that satisfy properties i) and ii), as the bosonic 
properties of the observed particles are not utilized. Such expansions are useful 
if one would like to study and characterize short range correlations in a general 
manner, and if no attempt is made to connect the structure of the shape of 
the BEC to the characteristics of the underlying space-time picture. (The 
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conditions for the convergence of such expansions are given in ref. [10].) 

The purpose of the present study is to provide a model-independent charac- 
terization of the BEC with a further additional assumption: 

Hi) The two-particle correlation function is related to a Fourier transformed 
space-time distribution of the source. 

This last assumption implies that the plane-wave approximation can be war- 
ranted, and additional possible short-range correlations due to e.g. the final 
state Coulomb or strong interactions as well as energy-momentum conserva- 
tion or the jet structure of the source can be removed from the data or their 
magnitude is under theoretical and/or experimental control. 

The purpose of the present paper is to determine in a general form the ap- 
proximate behaviour of the two-particle BEC under assumptions i) - Hi) and 
relate the parameters of the source functions to the parameters of the BEC in 
a generic form. 

1.1 Basics 

The two-particle BEC is defined as the ratio of the two-particle invariant mo- 
mentum distribution iV 2 (ki, k 2 ) to the product of the single-particle invariant 
momentum distributions iVi(kj) [11], where k stands for the three- momentum 
of the detected particle of mass m and energy E = y/m 2 + k 2 . The BECF is 
defined as 



for more details on the normalizations see refs. [5,11]. For simplicity, we assume 
that at large relative momenta, the two-particle correlation function tends to 1. 
Thus we neglect possible long-range correlations and normalization coefficients 
in eq. (1). 

As C 2 (ki,k 2 ) is a ratio of invariant momentum distributions, the correlation 
function is also an invariant observable. This property can be emphasized by 
using the C 2 (ki, k 2 ) notation, where k\ and k 2 are the momentum four- vectors 
of the particle pair. See also refs. [5,6] for more detailed discussions about the 
constraints between k\ and k 2 . 

We assume that C 2 (k±, k 2 ) is experimentally determined and satisfies proper- 
ties i) - Hi). These correlation functions can be evaluated at a certain fixed 
mean momentum K = 0.5(ki + k 2 ) as a function of various components of the 
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relative momentum q = ki — k2. The relative momentum can be decomposed 
into some experimentally preferred relative momentum components in 1, 2 or 
3 dimensions. For example, experimental data are presented in one dimen- 
sion in terms of the invariant momentum difference Qi = \J — (k\ — k 2 ) 2 , or 
as a function of various components of the relative momentum, for example, 
(?o, Qz, Qt) = (Ei - E 2 , k M - k 2 , z , yj(ki, x ~ k 2 , x ) 2 + (k hy - k 2<y ) 2 )■ 

Various directional decompositions of the relative momentum with respect to 
the mean momentum have been proposed. These decompositions include the 
Bertsch-Pratt, the Yano-Koonin-Podgoretskii and the invariant Buda-Lund 
decompositions. For further details, we refer the interested readers to the re- 
cent reviews on these topics, [5,6,12-14]. For our further considerations, the 
specific experimental choice for the directional decomposition of the relative 
momenta will not be essential. 

Essential background information on the generic structure of the peak of the 
BEC was summarized in ref. [15] regarding the impossibility of the Taylor 
expansion of the correlation functions around the Q = point, and in the 
appendix of ref. [16], where it was pointed out that the class of the stable 
distributions plays an important role in non-Gaussian BEC studies. The break- 
down of the Gaussian approximation around the Q = point was shown 
with the help of a theoretical analysis and a numerical simulation in ref. [17]. 
Explicit simulations, including long tailed and/or asymmetric source density 
distributions confirmed this result [18]. 

The above results are essentially the consequence of the fact that the two- 
particle phase-space volume grows as Q 3 for small values of the relative mo- 
mentum Q. Hence the number of pion pairs tends to zero as Q — > 0, which 
implies that the numerator and the denominator of eq. (1) tend to in the 
same limit, UmQ^ C 2 (Q) = "0/0". This is the reason why it is experimentally 
impossible to determine C(Q — 0) exactly: the experimental errors on the cor- 
relation function diverge as Q — > 0. (Not to mention and d j^p- 1 whose 
experimental errors diverge even more than that of C(Q)). So the two-particle 
correlation function can only be measured at non-zero values of Q > Q m in, at 
a scale set by the combination of factors such as the available statistics, the 
two-track resolution and the extent to which final-state interactions are under 
control. Typical values of Q m in are 5-10 MeV in present correlation measure- 
ments in heavy ion physics. The value of the Bose-Einstein correlation function 
at Q = can then be estimated with various extrapolation techniques, e.g. 
with the help of a (typically Gaussian) model assumption for the functional 
form of C(Q), or with the help of model-independent Edgeworth or Laguerre 
expansions [10]. 
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2 Bose-Einstein correlations for univariate stable distributions 



For clarity, let us consider first the case of a one-dimensional static boson 
source, and show the essential ideas in the simplest possible manner. (In the 
subsequent sections we generalize these results to more realistic, expanding 
sources that contain correlations between coordinate space and momentum 
space, and also to the case of more complicated derivations.) 

Thus we assume here that the source can be characterized by a factorized form 
of the emission function: 

S(x,k)=f(x)g(k), (2) 



where x and k are one-dimensional coordinate and momentum variables, re- 
spectively. The normalization is chosen such that 

fdxf(x) = l, jdkg(k) = (n), (3) 



where (n) stands for the mean multiplicity. The single-particle spectrum is 
obtained as 

N^k) = f dxS(x,k) = g{k). (4) 



Assuming chaotic particle emission, and the validity of the plane-wave ap- 
proximation, the two-particle Bose-Einstein symmetrized wave-function reads 
as 

ipk u k 2 (xi, x 2 ) = -7= [exp(ifcixi + ik 2 x 2 ) + exp(ikxx 2 + ik 2 xi)] . (5) 
a/2 



In the Yano-Koonin formalism [20], the two-particle momentum distribution 
is determined as 

N 2 (h,k 2 ) = J dx 1 dx 2 S(x 1 , h) S(x 2 ,k 2 ) |^ fcll fc 2 (xi,x 2 )| 2 , (6) 



and the resulting two-particle BEC is found to be 

C 2 (A; 1 ,A; 2 ) = l + |/(g 12 )| 2 , (7) 



which measures the absolute value squared Fourier transformed coordinate- 
space distibution function of the particle emitting source, where the Fourier 



5 



transformed source density (often referred to as the characteristic function) 
and the relative momentum are defined as 

/(912) = dx exp(iq 12 x) f(x), q 12 = h - fc 2 . (8) 



As the two-particle Bose-Einstein correlation function is determined by the 
modulus of this Fourier-transformed source distribution, but is insensitive to 
the phase of the Fourier transform, these modulus squared Fourier transforms 
measure only the relative coordinate distributions of the source, but they are 
not sensitive to shift parameters like the location of the center of the source 
distribution in space-time. It is frequently assumed (either explicitly or im- 
plicitly), that the function f(q) is analytic around q = and that its second 
order Taylor expansion characterizes its behaviour well even for large values 
of q. We shall see later that this imposes severe constraints on the possible 
choice of the source functions. As this point is essential, we emphasize that 
the analyticity of f(q) at q = is an assumption that greatly reduces the 
possible forms of the source densities. This may be clearly seen in the case of 
Levy stable distributions[21-23], which are characterized by four real param- 
eters, and thus form an uncountably infinite set. The most essential of these 
four parameters is the index of stability a, which may take any value in the 
interval (0,2]. However, the assumption that the f(q) characteristic function 
is analytic restricts this parameter to the special a = 2 case, corresponding to 
Gaussian Bose-Einstein correlation functions. 

If the function f(q) is an analytic function at q = and if the lhs of eq. (8) 
can be expanded into a convergent series, one obtains 

f(q)^l+iq(x)-q 2 (x 2 )/2 + ..., (9) 



where 



(x n ) = J dxx n f(x). (10) 



This, in turn, implies that 

C(q) = 1 + \]{q)\ 2 « 2 - q 2 ((x 2 ) - (x) 2 ) » 1 + exp(-g 2 J R 2 ), (11) 



where a Gaussian source radius parameter R is defined as the width of the 
source emission function, 

R = vV) - W- (12) 
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Such type of over-simplified derivations, repeated in multi-dimensional forms, 
are frequently utilized to "prove" that the Bose-Einstein correlation function 
"has to have a Gaussian form" , without investigating the domain of applica- 
bility of the approximations that yield this result. 

In fact, a Gaussian ansatz for the two-particle BEC frequently gives a good 
first approximation to the data, in particular, if the statitistical precision of 
the date is not adequate to perform a detailed shape analysis. One of the 
reasons for such a behaviour could be that the emission of elementary particles 
is a rather complicated, stochastic process, in particular in collisions with 
very high energies. If one assumes that there are many independent processes, 
that shift by 5x the coordinate x, and that the the final production point 
is a sum of many, similarly distributed, random shifts, x = J2i$ x i, an d the 
shifts are characterized by finite variances, then for such random variables the 
probability distribution for x tends to a Gaussian one, if all the conditions 
of the Central Limit Theorem are satisfied. As the Fourier transform of a 
Gaussian is also Gaussian in such cases the expected shape of the two-particle 
BEC is also a Gaussian. 

However, there are many random processes for additive random variables, 
where a limiting distribution exists, but the Gaussian version of the central 
limit theorem is not applicable. Such processes are characterized by large fluc- 
tuations, power-law like tails and the non-analytic behaviour of the character- 
istic function of the probability distribution for small values of its arguments. 
In mathematical statistics and probability theory, such distributions are fre- 
quently referred to as Levy, or stable distributions. 

Let us approach the problem of the shape of the two-particle BEC without a 
reference to a particular model, from the more general point of view of math- 
ematical statistics and the (generalized) Central Limit Theorems. Concrete, 
model specific applications shall be presented elsewhere [24]. 



3 Levy stable distributions and generalized central limit theorems 

In physics, as well as in the theory of probability, the probability distribution of 
a sum of a large number of random variables is one of the important problems, 
as such distributions are frequently realized in Nature. Limit distributions 
characterize the probability distributions of random processes in the limiting 
case when number of the elementary independent random subprocesses tends 
to infinity. 

Various forms of the Central Limit Theorem state, that under some condi- 
tions, the sum of a large number of random variables behaves as a Gaussian 
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distribution. Essentially, the conditions for the validity of the various central 
limit theorems are that the elementary distributions should have a finite mean 
and variance. The sum of a large number of such variables follows a Gaussian 
distribution, which in turn is a (frequently encountered) special case of limit 
distributions. 

For example, rapidity is a kinematic variable invariant under longitudinal 
Lorentz transformations. The final (observed) rapidity of a particle is a con- 
sequence of many random and independent actions, which can be considered 
as some rapidity shifts. In case of rescattering on a hadronic or on a partonic 
medium, the final rapidity of the particles is obviously a sum of many random 
rapidity shifts. If these subprocesses have finite means and variances, then 
the central limit theorem implies that the observable rapidity distribution has 
to be a Gaussian, regardless of the details of the probability distribution of 
rapidity shifts in the elementary collisions. 

The stable distributions are precisely those limit distributions that can oc- 
cur in Generalized Central Limit theorems. Their study was begun by the 
mathematician Paul Levy in the 1920's. A recent book by Zolotarev and 
Uchaikin [21] contains over 200 pages of applications of stable distributions 
in probabilistic models, correlated systems and fractals, anomalous diffusion 
and chaos, physics, radiophysics, astrophysics, stochastic algorithms, finan- 
cial applications, biology and geology. Stable distributions provide solutions 
to certain ordinary and fractional differential equations, and the breadht of 
their applications suggests that they can be considered as a class of special 
functions [21-23]. 



3. 1 Characteristic functions 

The Fourier transformed density distribution, eq. (8) is usually called the 
characteristic function in mathematical statistics. The stable distributions are 
frequently given in terms of their characteristic functions, as the explicit for- 
mulas of the corresponding source density distributions are known only in 
some special cases. 

The characteristic functions of univariate stable distributions can be defined 
using various convention schemes. In this work, we utilize the notation of 
Nolan [23]. In this work, the characteristic function of univariate stable dis- 
tributions is given in three different conventions, out of which we utilize the 
convention S(a, (3, 7, 5; 1) as this leads to natural generalization of data pa- 
rameterizations utilized before in Bose-Einstein correlations in high energy 
physics. In this convention, the characteristic function f(q), with other words, 
the Fourier-transform of a probability density given by the S(a, j3, 7, 5; 1) law, 
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has the following general form: 



f(q) = exp (- 7 Q |gr + t(3l a \q\ a sign(q) tan(— ) + iq5j , a ^ 1. (13) 

Thus univariate stable distributions are characterized by four parameters: the 
Levy index of stability (or characteristic exponent) < a < 2, the skewness 
parameter — 1 < /3 < 1, the scale parameter < 7 and the location parameter 
—00 < 5 < 00. We recommend to consult the book of J. Nolan for greater 
details and the a = 1 special case, [23], see in particular Fig. 1.3 of this work 
for an illustration of the stable densities in the S(a, (3, 7, 5; 1) parameteriza- 
tion [23]. 

Before discussing the univariate case in general, let us start with one of the 
simplest cases. Mathematically, this corresponds to the selection of the sym- 
metric stable distributions. Their relationship with the Gaussian case is also 
the most apparent. Let us also introduce physical units. Using the notation 
of Nolan [23], this choice corresponds to the S(a, (5 = 0, 7 = R/2~ , 5 = x ; 1) 
convention. 

For symmetric stable distributions, the Fourier transformed density distribu- 
tion has a simple form, 

f(q)=exp(iq6-\>yq\ a ), (14) 

where \z\ stands for the modulus (absolute value) of the variable z, and the 
support of the density function f(x) is (—00, 00). More complicated cases are 
discussed in the subsequent parts. In order to simplify the results, and to 
present results that are similar to the ones used in data fitting in high energy 
physics, it is useful to re-scale the scale parameter 7 of the Levy distributions 
and introduce a physical notation as follows: 

R = 2« 7 x = 5 (if /3 = 0)., (15) 

In the chosen S(a, j3, 7, 5; 1) notational system, for symmetric {(3 = 0) stable 
distributions, the parameter 5 coincides with xq, the location parameter of the 
distribution. 

A generic property of the stable distributions is that their characteristic func- 
tion is non-analytic at q = for all values of the index of stability a < 2. 
For symmetric distributions, the following small q expansion reflects this non- 
analytic behaviour: 

f(q)*l + iqx -±\qR\ a (16) 
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The analyticity of this expansion is restored only if a = 2, corresponding to 
the Gaussian limit. 

Levy-stable distributions yield the following, relatively simple form of the two- 
particle BEC: 

C(q;a) = l + ex P (-\qRn. (17) 

This form that has an additonal parameter, the index of stability a, as com- 
pared to the usual Gaussian (or exponential) distribution, where the value of 
a is fixed to 2 (or 1). 

Note that eq. (14) is not an approximation for the shape of the characteristic 
function around the q = point, but an exact and deep mathematical result, 
a consequence of the Generalized Central Limit Theorem, which gives the 
possible forms of the Fourier transform of the density distributions of the 
stable laws. Hence eq. (14) is valid for all values of q, and the corresponding 
form, eq. (17) for the two-particle BEC is well suited for a detailed shape 
analysis in the physically resolvable large q region. 

The simplicity of eq. (17) is the reason why we utilized the S(a, (3,7,6 — x ; 1) 
convention with a rescaling of the 7 parameter to R: In cases of a = 1 and 
a = 2 we obtain formulas that already had a number of applications in particle 
interferometry studies in high energy physics [5,6,14]. 

The physical meaning of the index of stability is that the symmetric stable 
distributions decrease for large values of the coordinate 1 as a power-law, 

f(x) oc x- 1 -*, for x > R, < a < 2, (18) 

while in the a = 2 special case corresponds to Gaussian source distributions. 

In general, stable distributions with < a < 2 are long tailed and are fre- 
quently related to the self-similarity of the generating mechanism. For exam- 
ple, in QCD, the theory of strong interactions, jets emerge within jets within 
jets etc, a process characterized by a self-similar branching and a fractal struc- 
ture, where the fractal dimension of the source is related to the anomalous 
dimension of QCD [26-28]. 

In the Lund model of hadronization, a self-similar string fragmentation process 
determines the distribution of the directly produced particles [29]. Among 
these, there is a big range of the distribution of the decay times of various 
resonances, and their decay may lead also to a power-law tail fluctuation of 
the effective source sizes [30]. Numerically, similar results were obtained by 
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Utyuzh, Wilk and Wlodarczyk in ref. [31] when considering Bose-Einstein 
correlations for sources with a fractal, power-law structure in space-time. 

A second order phase transition is expected if a heavy ion collision passes 
through the critical point of the quark-gluon plasma - hadronic matter phase 
transition, throught the end-point of the phase-boundary along which the 
transition is of first order. In the vincinity of the critical point the second order 
phase transition is expected to be signalled by power-law type of fluctuations 
in coordinate space. Coordinate distributions with power-law tails may obey 
the Generalized Central Limit Theorem, which essentially means that these 
density distributions, in the limit of large number of independent sources, tend 
to a Levy stable distribution with index of stability a < 2. 



3.1.1 Some examples 



1) The Gaussian distribution corresponds to the case of a 
parameter is R = 



2. The scale 

x 2 ) — (x) 2 , the standard deviation. The density distri- 
bution and the corresponding BEC are 



C(q) 



cxp 



(X 



(2nR 2 y/< 
1 + exp (-q 2 R 2 ) . 



2R 2 



— OO < X < oo, 



(19) 
(20) 



2) The Lorentzian (or Cauchy) distribution corresponds to the case of a = 
1. The scale parameter R is the standard scale parameter of the Lorentz 
distribution. The density distribution and the corresponding BEC are 



/(*) 
C(q) 



R 



71 R 2 + (x - x ) 2 ' 
1 + exp (-\q R\) . 



— OO < X < oo, 



(21) 
(22) 



Such a form was frequently applied to describe the Qi dependence of two- 
particle BEC in high energy heavy ion and particle physics [14,5]. 
3) The asymmetric Levy distribution is a slightly more complicated, but still 
relatively simple case, corresponding to a — 1/2,(3 = 1. The density has 
the following one-sided distribution 




/(*) = 



C(q) = l + exp(-J\qR 



8tt (x - x f/ 2 



cxp 



R 



8(x - x ) ) ' 



Xq < X < OO, 



(23) 
(24) 



This density distribution has been considered, although in a different for- 
malism, for the characterization of two-particle BEC in ref. [15]. 
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3.2 Asymmetric Levy distributions 



Adopting Nolan's S(a, f3, 7, 5; 1) convention for the univariate Levy stable 
densities with the rescaling of the scale parameter 7 to R, eq. (15) , the 
characteristic function f(q) has the following general form: 

/(?) = exp (~\R a \q\ a + ^iT|g| Q sign(g)tan(^) + iqx ^J , a + 1,(25) 

and we recommend to consult the book of J. Nolan for the a — 1 special 
case, [23]. 

This form of the Generalized Central Limit Theorem also implies that two- 
particle BEC-s, measuring the modulus square of the Fourier transformed 
source distribution, are insensitive to the asymmetry parameter f3 of the Levy 
stable distributions, as well as to the center x of the source. This is the di- 
rect consequence of the fact, that the two-particle Bose-Einstein correlations 
measure only the relative coordinate distributions, hence asymmetry or the 
position of the center of the distribution does not enter the observable corre- 
lations. 

If the source is characterized by a Levy stable distribution, the two-particle 
BEC can always be written in the form of eq. (17) in the whole range of q. 

Note that for (3 7^ 0, the interpretation of the 7 and 5 parameters of the 
S(a, (3, 7, 5; 1) is not straightforward, and this parameterization is not a con- 
tinuous function of the parameter a at the a = 1 point. See ref. [23] for a 
more detailed discussion on these points. 



3.3 Three-particle Bose-Einstein correlations for stable distributions 

Can one extract, at least in principle, the asymmetry parameter f3 of the stable 
laws, using the methods of particle interferometry? In order to answer this 
question, let us consider the three-particle Bose-Einstein correlation function. 
In general, it may have a complicated form [4,5]. However, if the particle 
emission is completely chaotic and if the plane-wave approximation can be 
warranted, then the three-particle BECF reads as follows: 

C 3 (l,2,3) = 1 + |/(1, 2)| 2 + |/(2, 3)| 2 + |/(3, 1)| 2 + 

+ 2^/(1, 2)/(2,3)/(3,l). (26) 
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where the symbolic notation f(i,j) = f(h — kj) = f{q%j) has been introduced 
to simplify the equation. 

The three-particle cumulant correlation function corresponds to the last term, 
k 3 (1, 2, 3) = 2ft/(l, 2)/(2, 3)/(3, 1), (27) 

where the two-particle cumulant correlation function is defined as 

« 2 (1,2) = |/(1,2)| 2 . (28) 

Hence, the normalized and symmetrized ratio 

"(1, 2, 3) = . "* (1 ' 2 ' 3) (29) 
2 > /k 2 (1,2)k 2 (2,3)k 2 (3,1) 



turns out to be a simple function of the asymmetry parameter of the Levy 
stable distributions: 



w(l,2,3)=cos |iTtan(— )[£ Irfsign^U for a ^ 1, (30) 

(for the special case of a = 1, see again ref. [23].) In the above equation, 
the summation is taken over the cyclic permutations, = (1,2), (2,3), or 
(3, 1). Note that the displacement parameter 5 = x cancels from this result, 
too. However, this parameter, related to the center of the particle emitting 
source, is not considered as an essential quantity in the physics of particle 
interferometry. For simplicity we may assume that the origin of our coordinate 
system is fixed by the condition xo = 0. 

Let us emphasize that the result of eq. (30) is valid only within the plane-wave 
approximation and neglecting possible partial coherence in the source [4,5]. If 
these conditions hold, the asymmetry parameter — 1 < /3 < 1 can be de- 
termined from the relative momentum dependence of the normalized three- 
particle cumulant correlation function w, which is denoted also [32,33] as 

w = cos(0). (31) 



The angle is directly proportional to (3, the asymmetry parameter of the 
stable distributions, and the coefficents of proportionality can be determined 
from the two-particle correlation data. Thus we have proven, that for a one- 
dimensional Levy stable distribution the asymmetry parameter of the source 
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can in principle be determined from a combined use of the two- and three- 
particle BEC-s. Thus all the essential parameters of the source density can be 
reconstructed within this class of density distributions. 

We would like to note, that the cases of f3 — +1 and (3 = —1 are particularly 
interesting, as in these cases the support of f(x) is one-sided, corresponding 
to [5, oo) and (— oo,<5], respectively. Such kinds of Levy stable distribitions 
may have applications e.g. in the description of the temporal processes in the 
collisions of particles at high energies: particle production cannot start before 
the initial collision - hence the corresponding time interval has to be limited 
from below. 

After its parameters a,j3 and 7 are determined, a Levy stable density distri- 
bution f(x) can be straightforwardly evaluated using public domain numerical 
packages [36]. In most of the cases these functions are also known in closed 
form. The difficulty is that they can be expressed in terms of very general 
kind of functions, Fox's H-functions, see Schneider's papers for greater de- 
tails [37,38]. 

In this sub-section we neglected a possible core-halo correction, because it is 
straightforward to shown that it cancels from the final result on w, eq. (30). 

Note, that in the more general case of expanding sources, Bose-Einstein cor- 
relations are sensitive only to the so called lengths of homogeneity, which 
correspond to the local density distribution for particles that have a given 
momentum k. In order to reconstruct the complete picture of the density dis- 
tribution for all particles, one has to perform a combined analysis of the single 
particle spectra and the two-particle BEC-s [16,39,5]. 

For such type of expanding sources, the radius parameters are known to de- 
crease with increasing (transverse) mass of the particles [5,6]. This implies, 
that the sensitivity of three-particle correlations for the shape parameters of 
the source distributions is enhanced when the heavier kaons are utilized in- 
stead of the pions, as follows if one evaluates eq. (30) using a reduced radius 
parameter. 



4 Multivariate symmetric Levy distributions 

Multivariate Levy stable distributions have properties similar to the uni- 
variate case, their study is even presently a research topic in the theory of 
mathematical statistics, see ref . [40] . In this section, we assume that the source 
is symmetric, to simplify the presentation. 
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Let us denote the coordinate 3-vector by x, and the time by the variable t, 
while a space-time four- vector is x = (t, x) , a momentum space four- vector is 
k = (E,k). 



The simplest case corresponds to a static (non-expanding) source, that emits 
all the particles instantaneously. This is described by a factorized form of the 
emission function, 



//the density distribution /(x) is a multivariate stable distribution that fea- 
tures ellipsoidal contour lines, than the reconstruction of the source density 
distribution becomes a one-dimensional, solvable problem [25] . 

The corresponding three-dimensional, symmetric, sub-Gaussian Levy-stable 
characteristic function is 



where x stands for the center of particle emission, and i?? ■ corresponds to a 
symmetric matrix R it j = Rj ti with dimensions of radius squared; a is the index 
of stability similarly to the one-dimensional case, and q = ki — k 2 stands for 
the relative momentum. 

Again, we have applied a rescaling as compared to the S(a, 0, 7, 0; 1) class 
of sub-Gaussian distributions, that are discussed in ref. [25]. The rescaling 
corresponds to the application of eq. (15) interpreting that equation in a ma- 
trix notation. We chose a notation that is adopted to the convention used in 
correlation studies in high energy physics, so that the Gaussian limiting case, 
a — 2, corresponds to the 



conventional multivariate Gaussian forms. 

Hence, for symmetric, sub-Gaussian multivariate Levy distributions, the two- 
particle BEC has the following simple form 



S(x,k)=f(x)g(k)6(t-t ). 



(32) 




(33) 




(34) 




(35) 
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The components of the squared radius matrix are symmetric, R 2 j = i? 2 ,. The 
components of this matrix and that of the relative momentum can be deter- 
mined in the experimentalists' favourite reference frame. The point is that 
QiR'ijQj is a (non-relativistic) bilinear form, which is (non-relativistically) 
invariant for the choice of the coordinate system, hence qR 2 q can be decom- 
posed in any frame. 

The ellipsoids are characterized by the inverse of the R 2 matrix, which is 
denoted symbolically by R~ 2 , and frequently called the covariance matrix: 

£4^ = ^. (36) 

j 



The large distance behaviour of the density distribution can be analytically 
determined as follows: 



/(x)^(^i?- 2 ^,)-^, (37) 



so the source density decreases as a power-law for sufficently large values of 
the arguments, with < a < 2. (The case of a = 2 corresponds to Gaussian 
source distributions). 

The elliptically contoured density distributions can be calculated by integrat- 
ing a one-dimensional distribution [25]. The ellipsoidal contour lines can be 
transformed to spherical contours by realizing that the two-particle BECF is 
characterized by R 2 , a symmetric matrix {R 2 )ij = {R 2 )ji , which always can 
be written as R 2 = R-R T . The matrix R can be considered as a sort of 'square 
root' of the squared radius matrix. Physically, the matrix R gives the orienta- 
tion of the prinicipal axis of the ellipsoids of the density profile in coordinate 
space, and characterizes also the stretchings along the principal axis that are 
needed to transform a unit sphere to an ellipsoid corresponding to the contour 
lines. Hence the bilinear form in the exponent of the BECF can also be written 
as qR 2 q = |Rq| 2 . 

The density distribution is expressed, in d > 1 dimensions, for a < 2, by eqs. 
(1) and (6) of ref. [34], 



oo 

/a(s(x)) = (27rdetiff)'/* / md ~ l (ts(x)) "" /2 J'/t-iM*))*-*', (38) 
s(x) = |i2 _1 x| (39) 

where Ji(z) is the i^ 1 order Bessel function of the first kind, the matrix R^ 1 
is the inverse of the radius matrix R, and the function s(x) corresponds to 
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a scaling variable that is constant over the ellipsoidal surfaces characterizing 
the constant source densities. It is relatively easy to see the physical meaning 
of this equation in d = 3 dimensions, in the system of ellipsoidal expansions, 
where the radius matrix corresponds to it! = diag(X, Y, Z), its inverse is R^ 1 = 
diag(l/X, 1/Y, l/Z) and the scaling variable of the ellipsoidally symmetric 
source distribution is given by 



Note that such scaling variables arise naturally in ellipsoidally symmetric solu- 
tions of fireball hydrodynamics, that may correspond to non-central collisions 
of heavy ions in the non-relativistic energy domain [41-43]. 

In addition to the above analytic result, eq. (38), the case of multivariate sym- 
metric Levy laws is discussed in great details in refs. [34] and [35]. [We don't 
recapitulate here the well-known results for multivariate Fourier-transforms 
of Gaussian densities, corresponding to the case of a — 2.] From the point of 
view of the applications, it is important that public domain software packages 
are available [36] , that can be utilized to plot the density distribution of stable 
laws as a function of the scaling variable s, if the parameter a is determined 
from a multi-dimensional fit to the two-particle BEC-s. 

Multi-dimensional BEC-s are frequently studied in high energy particle and 
heavy ion physics. For the sake of completeness, and the applicability of our 
results to these cases, we give the two-particle BECF in a general form, using 
the Wigner- function formalism of refs. [44,13,5]. 



4-1 Wigner-function formalism 

From now on, x and k stand for four-vectors. Using this notation, most of the 
results for the univariate case can be generalized to the relativistic case. 

In the Wigner-function formalism, the particle emission is characterized by the 
source or emission function S(x, k). In the non-relativistic limit, this function 
corresponds to the time- derivative of the Wigner function, which itself is the 
quantum mechanical analogy to the classical phase-space distribution and this 
S(x, k) can also be considered as a relativistic, covariant form of the Wigner 
function. 

The invariant momentum distribution and the two-particle BEC can be ex- 
pressed in the Wigner-function formalism as follows: 




(40) 
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N 1 (k 1 ) = J d 4 xS(x,h), (41) 

C 2 ( kl , fe) = 1 + . « 1 + S^, (42) 

21 2J 5(0,^)5(0,^) |5(0,K)|2' 1 ; 

X) = /" d 4 a; exp(^)S(x, K), (43) 

^ = |(fci + A; 2 ), q = k 1 -k 2 . (44) 

By introducing the notation 

= WT^-y < 45 > 



the results of the earlier section are recovered, but with all the parameters of 
the effective source depending on the mean momentum K. 



4-1.1 Core-halo correction 

If some of the particles are emitted from the decays of very long lived reso- 
nances, such that their decay width is much smaller than the typical 5-10 MeV 
width of the first bin of the experimentally determined two-particle BECF, 
then the decay products of these particles contribute to a halo that cannot 
be resolved with the help of two-particle BEC measurements [1]. Within the 
core-halo picture, their effects can be taken into account by assuming that the 
single-particle emission function can be written as a sum of a resolvable core 
contribution and an unresolvable halo contribution, 

S(x,k) = S c (x,k) + S h (x,k) (46) 



and as a result, an intercept parameter A is introduced to the two-particle 
BEC, 

C 2 (h, k 2 )^l + \ \ S A q,K ^\ (47) 
1 ; \S c (0,K)\^ 1 ; 



where the effective intercept parameter < A = A(k) < 1 is interpreted as 
a (mean-momentum dependent) measure of the fraction of bosons emitted 
directly from the core: 

A _ r n c (k) 

[N c (K) + N h (K) 



(48) 
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and the emission function of the core determines the relative-momentum de- 
pendence of the BECF. 



This formalism is relevant in case of correlation studies for expanding systems. 
Assuming that the local regions of homogeneity can be characterized by a 
symmetric multivariate Levy stable distribution, the general form for the two- 
particle BECF is as follows: 



C 2 (h, k 2 ) — 1 + Aexp 



3 

E R 



ijQiQj 



.a/2 



(49) 



where all the parameters may depend on the mean momentum of the pair, 
(A,i%,a) = (A(k),i%(K),a(K)). 

It is at this point where the large q behaviour of the Fourier transformed core 
density distribution becomes essential: within the core-halo picture, the ex- 
act intercept parameter of the two-particle BECF is A xc t = C 2 (k, k) — 1 = 1, 
however, it is assumed that the halo of long lived resonances creates an unre- 
solvable peak in the q < g m i n region, where the two-particle relative momentum 
resolution is typically g min ~ 5 MeV/c. This peak makes the analytic expan- 
sions around q = impossible, and the large q behaviour of the core correlator 
becomes important. Fortunately, this is uniquely determined with the help of 
the index of stability a for the chosen class of stable distributions. 

The resulting BEC-s all have a generalized exponential form, but, as com- 
pared to the usual Gaussian case, there exist an additional scale transforma- 
tion qi — > q^ 2 which sometimes is also referred to as stretching. Hence, the 
resulting correlation functions are multivariate stretched exponential param- 
eterizations. For a = 2, the earlier, Gaussian forms are recovered, while for 
a < 2 the correlation function becomes more peaked than a Gaussian, and it 
develops longer tails as well. As such a peaked, but almost Gaussian behaviour 
is frequently seen in the experimentally determined Bose-Einstein correlation 
data in high energy particle and heavy ion induced reactions, we may hope 
that the above stretched exponential or sub-Gaussian parameterization pro- 
vide a useful tool for the characterization of such datasets even in the case of 
the multi-dimensional analysis. 

It follows from the above discussion, that perhaps the presentation of the two- 
particle Bose-Einstein correlation data in the usual, linear-linear plot of C(q) 
versus q is not particularly useful. In order to exhibit the possible power-law 
structure of the exponent of the correlator, a log-log plot of the two-particle 
cumulant correlation function, In {ln[C(g) — 1] — In A} versus In q would be the 
most appropriate. On such plots, the two-particle correlation functions would 
appear as straight lines and the slope parameter of the lines would be given 
by the index of stability —a. 
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Below let us give two examples, where we highlight the essential model as- 
sumptions and show the result. For more detailed discussion on these topics, 
see refs. [5,6]. 



4-2 Levy stable distributions for non-relativistic, expanding sources 



Let us consider a popular, but not boost-invariant parameterization of the 
two-particle BEC-s, the Bertsch-Pratt [45,46] or side, out, long decomposition, 
where q = {q s ,q ,qi), where qi = ki >z — k<i, z is the longitudinal, q Q = (q x K x + 
q y K y )j ^Kl + is the out (ward) and q s = (q x K y - q y K x )/ + K% is the 
side(ward) relative momentum komponent. This decomposition is frequently 
studied in the LCMS frame [47], where K z = 0, because in this frame, the 
temporal component of the source couples only to the out direction, see e.g. [5] 
for details. This can be shown by utilizing the mass-shell constraints for the 
four-momenta, from which it follows that in the LCMS 

Hi l + Ej2 



If correlations between coordinate and time are negligible, a model assumption 
may be that the emission function at every value of the particle momenta 
can be written as a factorized coordinate-space distribution and a temporal 
distribution, 

S(x, k) = S(x, k)h(t) (51) 



and we may assume that at every value of k the coordinate space distribu- 
tion corresponds to a three-dimensional, sub-Gaussian Levy distribution with 
index of stability a x , while the temporal distribution would correspond to a 
univariate Levy distribution with index of stability at, and in general these 
two indexes may be different from each other. 

Such kind of factorization assumptions can be approximately made not only for 
non-relativistic expanding sources but sometimes also for relativistic expand- 
ing sources when the emission function is evaluated in the LCMS frame [47], 
and the motion is non-relativistic with respect to this frame. 

The resulting two-particle BECF for cylindrically symmetric sources can be 
written as 



C 2 (ki,k 2 ) = 1 + Aexp 



i,j=s,o,l 



(52) 
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An important feature of eq. (52) is that the spatial and the temporal scales may 
enter the relative momentum dependence in the out direction with different 
powers, a x and a t , hence, a detailed shape analysis of the BECF in terms 
of BP variables can be utilized to extract these different indexes as well as 
the information on both the spatial and the temporal scales. However, as is 
well known, the at = a x case becomes degenerate, and only the combination 
R 2 + fl 2 At 2 can be determined from the data analysis in this exactly Gaussian 
case. 

A more realistic case is to include possible correlations between production 
time and position, which goes beyond the above factorized ansatz. 



4-3 Levy stable distributions for the invariant Buda-Lund variables 



For collisions with very high energy, the particle emission process becames a 
highly relativistic phenomena. In this case, the invariance of the emission func- 
tion can be reflected if the longitudinally boost-invariant proper-time variable 
r = \Jt 2 — r 2 z is utilized, and the space-time rapidity rj = 0.5 log[(t+r 2 ) / (t—r z )] 
is also introduced as a hyperbolic, boost-additive coordinate. 

In a factorized form the Buda-Lund (BL) parameterization assumes the fol- 
lowing structure for the emission function [16]: 

S(x,k) = H*(T)G*(r))I*(r x ,r y ), (53) 



where the subscript * denotes an implicit momentum dependence. The effec- 
tive proper-time and space-time rapidity distributions H*(r) and G*{r)) are 
assumed to have a uni-variate Levy distributions with indexes of stability a= 
and a\\, while i*(r x , r y ) may have a bivariate Levy distribution with index a±. 
Here we refer to the symbolic notation of the invariant temporal and the par- 
allel relative momentum differences, Q= and Q\\, being conjugated variables 
to the space-time variables (r, rj). 

The corresponding factorized Levy-stable BL parameterization of the two- 
particle BECF reads as 



C7(fci,A; 2 ) = l + Aexp(-| J R = g = | a = - \R\\Q\p - \R±Q±\ a± ) , (54) 

Q= = m t) i cosh(yi - rj)- m t>2 cosh(y 2 - rj) (55) 

Q\\ = mt,i sinh(yi - rj) - m tj2 sinh(y 2 - rj) (56) 

Q± = y/Qi + Ql- (57) 
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In these equations, m t ,i = ym 2 + kf is the transverse mass, yi = 0.51n[(_Ej + 
k Z;i )/(E i — k Z;i )] is the rapidity of particle i and TJ stands for the space-time 
rapidity of the point of maximum emissivity for particles with a given fixed 
four- momentum fcj. Note that the above equations for the definition of the 
temporal and parallel invariant momentum differences are equivalent to the 
ones given in ref. [5], but the invariant decomposition is perhaps more apparent 
in the above form. 

The three different indexes of stability, a=, a\\ and a± appear as new fit 
parameters, due to the assumed factorization property of the boost-invariant 
temporal, longitudinal and transverse distributions from each other. All are 
Levy indexes for the corresponding spatial distributions, hence satisfy the 
usual unequality < CKj < 2 for all j = (=, \\, \\)- 

All the five scale parameters, A, R=, R\\, R± and rj may depend on the value of 
the mean momentum K, due to the expanding nature of the particle emitting 
source. 



5 Applications to NA22 and UA1 data 

The case of the uni-variate stable distributions can be checked, relatively easily, 
against the precision data on the two-particle short-range correlation functions 
as measured in the CERN experiments NA22 [48] and UA1 [49]. 

We emphasize, that we are not attempting here to reconstruct the detailed 
space-time structure of the particle emission process, or, interpret the result 
in terms of a physical process that may result a Levy stable law. We merely 
consider this application as a first test of the feasibility of our method as a 
quantitative tool for the analysis of two-particle BEC in high energy physics. 

The parameters for the best Levy stable fits are summarized in Table 1. 

The two-particle BEC for Levy stable distributions is a stretched exponential. 
Hence the title of Figure 1 where the best Levy stable fits are shown for both 
the NA22 and the UAl data sets. 

The plots of the best Levy fits are clearly better than even the best exponential 
fits to the same curve (not mentioning the "best Gaussian" approximations, 
who miss completely the strongly rising peak near to the Q = point). 

The quality of the fits is approximately similar to those obtained with the 
help of Laguerre expansions. However, an advantage of the Levy stable laws 
and fits is that they relate the exponent a to the power-law tail behaviour of 
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density fluctuations in the large coordinate-space distributions of elementary 
particles. We discuss a possible interpretation of our result in terms of QCD 
cascades in a separate paper [24]. 

The following statements can be made without any reference to a dynamical 
model: 

Both data sets can be described by the Levy stable fits in a statistically ac- 
ceptable manner. In case of the NA22 data, the confidence level of the fit is 
about 80 % even when one of the fit parameters is fixed, as determined from 
a Laguerre expansion method in ref. [10]. Without fixing this parameter, the 
minimum was found to be degenerate and the confidence level even higher. In 
case of the UAl data, the confidence level of the fit is only 6 %, a rather low 
value. 

In case of the NA22 data set, the fitted value of the A parameter is within 
errors close to the maximum value of 1 which is a theoretical upper limit 
within the core-halo picture. Such an upper limit is violated by the fit to the 
UAl data set. This may indicate that in the UAl data there some non-Bose- 
Einstein type of short range correlations present. However, fixing A = 1 a 
Levy stable fit to the UAl data becomes still possible, but only with a small a 
confidence level of about 1 %. Fixing the normalization to that of the Laguerre 
expansion for the same function, ref. [10], a Levy stable fit becomes possible 
with only a very small, 0.1 % confidence level. In all of the above cases, a 
remains significantly below 1. 

Fixing a — 1 (exponential fit), the confidence level of the fit to the NA22 
data decreases from 80% to 24 %, in case of the UAl data, from 6 % to 10 -8 . 
Fixing a = 2 (Gaussian fit), the confidence level of the best fit to the NA22 
data decreases to 10~ 8 , while it becomes essentially in case of the UAl data 
set. 

The conclusion is that Levy stable (or stretched exponential) fits describe the 
NA22 and the UAl data on short-range correlations of particles in high energy 
physics. In case of the NA22 data, the confidence level of the fit is very high 
and the fitted parameters can be interpreted within the core-halo picture, with 
a fully resolved halo. In case of the UAl data, the best fits are obtained when 
the parameters were utilized outside of their usual domain. This result may 
indicate the presence of short-range correlations of not Bose-Einstein type in 
the UAl data set. 



23 



6 Discussion and conclusion 



Levy stable laws appear in many branches of physics, mathematics, biology, 
economy, computer science and other areas, where the scale of fluctuations may 
be characterized by long tails and an asymptotic power-law like behaviour. As 
the Levy stable distributions satisfy a generalized central limit theorem, they 
may appear naturally also in case of high energy elementary particle physics, 
where self-similar branching processes are known to exist on the large momen- 
tum transfer scale, a self-similar string fragmentation process is expected to 
describe the distributions of directly produced hadrons, and where the distri- 
bution of abundances of various resonances may also follow an approximate 
power-law like behaviour as a function of the life-time of the resonance. Similar 
processes are expected to exist also in high energy heavy ion collisions, where 
additional fluctuations from a hadronic or quark matter may also modify the 
spectrum of these fluctuations. 

We have argued, that additive stochastic processes may determine the space- 
time density distribution of elementary production in highly energetic reac- 
tions of particles or heavy ions. In this case, if a generalized Central Limit 
Theorem is applicable, the probability distribution of particle emission corre- 
sponds to a(n univariate or multivariate) Levy stable distribution. 

Although the characterization of the Levy stable distributions is not straight- 
forward, by now these distributions are numerically as well as analytically 
accessible. Their most important property is that the Fourier transformed 
emission function has in general a non-analytic structure in the vanishing rel- 
ative momentum region. 

We have shown, that due to this reason, the two-particle Bose-Einstein correla- 
tion functions are characterized by a sub-Gaussian two-particle Bose-Einstein 
correlation function, with an index of stability < a < 2. The analytic be- 
haviour and the Gaussian limiting case corresponds to the a = 2 special 
case only, a single point out of the uncountably infinite possibilities in the 
parameter space. Hence there is no reason to believe that the two-particle 
Bose-Einstein correlation functions have to behave as (multivariate) Gaussian 
forms. 

We have demonstrated, in the simplest one-dimensional case, that the NA22 
and the UAl two-particle Bose-Einstein correlation functions are well de- 
scribed with a Levy-stable law and the index of stability was found to be 
significantly smaller than 2, i.e. the deviation from a Gaussian behaviour in 
this case can be quantified and characterized with our method. 

We have also determined the form of the two-particle Bose-Einstein corre- 
lation function in case the particle emission is characterized by a multivari- 
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ate Levy-stable law. In particular, we have generalized the non-relativistic 
Bertsch-Pratt and the relativistically invariant Buda-Lund parameterizations 
for multivariate Levy-stable laws. 
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Table 1 



Best fits to UA1 and NA22 two-particle correlations using a Levy stable law 





UA1 


NA22 


Parameter 


Value 


Error 


Value 


Error 


M 


1.34 


± 0.03 


0.95 (fixed) 




A 


1.85 


± 0.07 


1.15 


± 0.17 


R [fm] 


6.36 


± 0.33 


1.33 


± 0.30 


a 


0.49 


± 0.01 


0.67 


± 0.07 


X 2 /NDF 


56.5/42 = 


1.35 


27.67/35 = 0.769 




CL 


6.6 % 




80.6% 
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Stretched exponential fit 




Fig. 1. The figures show -D| which is proportional to the two-particle short range cor- 
relation function, as measured by the UAl and the NA22 collaborations at CERN. 
The lines stand for the stretched exponential (or Levy) fits, which are able to re- 
produce the data with an acceptable x 2 /NDF, while the dashed lines stand for the 
result of the best exponential fits. The fit results are summarized in Table 1. 
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